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A device model of biological molecule sensors based on semiconductor nanowires has been developed. This 
model of a bioFET is based on the concept of the electrolytic absolute electrode potential. From that starting 
point a semiconductor device model of the nanowire solution biomolecule system was derived. The model 
includes the Gouy-Chapman-Stem model of the salt solution double layer, site binding charges on the electrode 
surface, and biological molecules in the form of a membrane layer. A simple method of solving this model 
is presented using the finite element method. Some examples showing the general properties of the model are 
given. 



I. INTRODUCTION 

The idea of using a field effect transistor (FET) to de- 
tect charged molecules in a solution was first introduced by 
Bergveld. ^ ^ Since then there has been much research into us- 
ing FET based devices to detect charged biological molecules 
(bioFET).^'* In the past 10 years or so research has turned 
to nanowire FETs, with a number of workers reporting high 
sensitivities,^ "*^ but the results are still variable and not well 
understood."* 

In order to understand the experimental results, there is a 
need for a simple, accurate, and easy to use model. The basic 
principle of bioFET operation is that charged ions attached to 
the gate oxide will attract or repel carriers in the FET channel, 
changing the channel conductivity. A model needs to describe 
the structure of the mobile ions in the electrolyte, surface 
charging of the gate insulator, along with the biomolecules 
and how these interact with the carriers in the semiconductor. 

A number of analytic treatments have been published.^El 
These use approximations to render the mathematics solvable, 
and are necessarily incomplete, so only apply to restricted 
cases. The complexity of the problem leads to using numer- 
ical methods. Early numerical models were one dimensional 
and divided the system into a number of layers; essentially a 
series of capacitance's. Poisson's equation was solved in each 
layer in sequence while matching boundary conditions. This 
process was then i terated until an over all consistent solution 
was obtained.nSGI] This is a somewhat awkward and complex 
procedure that has only been applied to planar structures. 

Semiconductor style modeling is the next step in numeri- 
cal modeling.nSEl This merging of semiconductor and elec- 
trolyte regions is a complex system, leading to examples of 
incompletely explained or even incoiTect models. Recently, 
Dutton published results of a fairly complete numerical de- 
vice model. Even this publication does not explain the basis 
of the model, boundary conditions, or how they performed 
the calculations. Without these details, it is difficult for other 
workers to try to duplicate their results. A very detailed model 
solves the full semiconductor model in the Si, including the 
continuity equations and a separate Monte Carlo simulation 
of the charge in the layer of biomolecules.'^ This model is 
complicated and the calculations are difficult. Since the semi- 
conductor part and the biomolecule part are separated, the en- 
tire solution may not be self consistant. 



Other recent worlP uses proprietary software to solve the 
semiconductor equations in 3 dimensions. This allows sim- 
ulation of incomplete coverage of a bioFET using charged 
cubes to represent biomolecules. This model, however, does 
not fully represent a bioFET as it does not include features 
such as the Stern layer and specific surface charges (eg. site- 
binding model). 

In many experimental situations it can be arranged that the 
bioFET surface is uniformly coated with biomolecules. This 
case can be modeled as a 2 dimensional membrane, which will 
be done in this paper It would seem expedient to concentrate 
on this simpler situation, at least until basic bioFET behaviour 
is understood. Then one could move to the more complicated 
case of partial coverage. 

The purpose of this paper is to give a basic bioFET de- 
vice model of a semiconductor-interface-electrolyte system 
(Si/Si02/Sol) and a simple method to calculate the results. 
This model should be based on semiconductor physics and 
electrochemical principles so that it can serve as a basis for 
further development. The idea is to include all the important 
features of the bioFET system and using as simple a calcu- 
lation as possible. This is to allow easy comparison between 
theory and experiment to facilitate investation of the proper- 
ties of bioFET devices. 

The semiconductor-interface part will be represented by Si 
with a thin Si02 surface layer and the electrolyte as a salt solu- 
tion. The next section will discuss the electrolytic cell in terms 
of absolute electrode potentials. In the following section the 
view will change from the electrochemical view to the semi- 
conductor device view. The device model will include a band 
structure for the Si/Si02/Sol system, electron and hole den- 
sities in the Si and the structure of ion concentrations in the 
solution. A membrane model will be used to model a charged 
biomolecule layer. 

A finite element method of calculation will be outlined. 
The results for some basic examples will be used to demon- 
strate cylindrical nanowire end conditions, and calculation of 
response or sensitivity. 

II. ELECTROCHEMICAL CELL 

While this paper is concerned with modeling the 
Si/Si02/Sol system, an experimental system must consist of 
a complete electrolytic cell with two electrodes. The concept 
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of absolute electrode potential will be used here following the 
detailed work of Trasatti and Reiss.'** Consider the cell 



1 2 3 
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(1) 



The Si I Sol interface represents the system that we wish to 
model. The thin Si02 layer is omitted here for simplicity. The 
Sol I M junction represents a reference electrode. A complete 
circuit is formed by connecting the reference electrode, M, to 
an ohmic contact on the Si electrode with a wire made of the 
same metal as the reference electrode, M'. An ohmic contact 
means that M' and Si are in electronic equilibrium, thus the 
electrochemical potentials in M' and the Si are equal. The full 
cell potential can be written as the difference of the absolute 
electrode potentials, 



E - E^l - E^ 



(2) 



There are a number of possible definitions of absolute elec- 
trode potential. The definition that is useful here uses a "free" 
electron in the solution as the reference state, Trasatti's \E^^ 
For the Si electrode, this gives 
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where /Ig' and /I^"' are the electrochemical potentials of elec- 
trons in the Si and the solution, respectively. Substituting the 
definition of electrochemical potential gives 
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where and IX^°^ are the chemical potentials of electrons in 
the Si and the solution, respectively, is the inner electro- 
static potential in the Si bulk and 0^°' is the inner electrostatic 
potential in the solution bulk, q is the elemental charge. The 
inner potential is given ^ ~ X^W^ where % is the surface 
polarization and i// is the outer electrostatic potential, is 
generated by free charges at the surface of a phase plus ex- 
ternal field sources. The electron work function of a material, 
<I>, is the negative of the electron real potential, a^, given by 
= —as. = —jXe+qX- Now Q can be written in terms of 
measurable quantities. 



pSi _ 
^abs ^ 



.Sol 



q 



q 



^Si_^Sol 



(5) 



In the case of electronic equilibrium for the Si Sol electrode. 



£|s = Oand 



Ay/ = v/^' - 



^Si (j)Sol 

q q 



(6) 



This gives the boundary condition for the electrostatic po- 
tential. For a detailed picture of the inner potential in the 
Si/Si02/Sol system, including the polarization at the inter- 
faces see the paper by Bousse.*^ 



The voltage on the Si with respect to the solution, V = £. 



can be expressed using (|2]) as V = E^^ + 



i-app: 



abs' 

where Vapp is 



applied between the electrodes by an external power supply. 

The approach to this problem used by Ref. [16] is to ignore 
the chemical potential of the reference electrode, ie. to as- 
sume an ohmic contact to the solution. Then, they moldel the 
solution as a semiconductor with a 1. 1 eV energy gap and with 
an artificially defined electron affinity. The value of which is 
determined by matching their model calculations to experi- 
mental data. 



III. DEVICE MODEL 

The model of the system Si/Si02/Sol will start by con- 
sidering the equilibrium state and then be extended to non- 
equilibrium with essentially zero current. It is important to 
note that, it is common, in semiconductor modeling, to use 
the outer potential, y/, rather than the inner potential, 0. y/ 
does not include the surface and interface polarization and so 
is continuous at interfaces. The effects of surface polarization, 
X, is, arbitrarily, included in the value of the work functions. 

In the semiconductor device model one must solve the Pois- 



son equation,^ 



I22I23I 



(7) 



e is the position dependent dielectric constant and is the 
electric potential, p is the charge density and is composed of 
electron, e, hole, h, and dopant densities in the Si and of ions 
in the solution. 

In the general, non-equilibrium case, continuity equations 
for each mobile particle must also be solved. In the model 
presented here a voltage applied to the Si with respect to the 
solution will result in essentially zero current due to the insu- 
lating Si02 layer. With zero current, the particle conservation 
equations are trivial and the applied voltage can be accounted 
for in the expressions for the charge densities. 

In the next section the band structure of the system will be 
described, followed by a development of the charge carrier 
statistics. Then the structure in the electrolyte will be pre- 
sented. 



A. Band structure 

The Si/Si02/Sol system will be treated as a semiconduc- 
tor he teroj unction system.^^l The solution can be thought of 
as a low carrier density metal or a very small band gap semi- 
conductor Consider isolated p-type Si and solution phases 
as drawn in Fig. [T^). The Fermi level in the solution coin- 
cides with the energy of an electron in the neutral solution. 
The work function for each material is the energy difference 
between the Fermi level and the vacuum level. 

Fig.[T]3) shows a band diagram where the system is at elec- 
tronic equilibrium so that the Fermi level is constant across 
the system. To achieve this, charge must be redistributed such 
that the resulting electrostatic potential, \//, bends the vacuum 
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FIG. 1: Energy diagrams, a) Isolated p-type Si and solution phases 
showing the work functions, "t^' and as well as the Fermi lev- 
els, Ef , and intrinsic levels, E,-. In the solution, the intrinsic level 
is equal to the Fermi energy, b) Si-Si02 -electrolyte system at equi- 
librium. The charge redistribution generates an electric field which 
changes the energy of the vacuum level, Evac, and the bands. 



energy level by an amount equal to the difference of the work 
functions. This is the builtin potential and is expressed by (|6|. 

When a voltage, V, is applied to the Si, the potential, y/, is 
increased by V and the electron energy bands and the Fermi 
level are lowered by —V . The difference in the bulk values of 
y/ in (|6]l is modified to 



4)Si 



(j)Sol 
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(8) 



The work function of neutral Si depends on the dopant density, 
and the intrinsic energy. It can be expressed \>y^ 



-kT\n 



E/r = E, — kTln 
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n — type (9a) 
p — type (9b) 



where k is Boltzmann's constant, T is the absolute tempera- 
ture, n, is the Si intrinsic carrier density, Nd is the donor con- 
centration, and Na is the acceptor concentration. A parameter 
can be defined, which is the difference between the Si intrinsic 
level and the electrolyte neutral level. 



5,=Ef-E|°'=E, + <I>^°i. 
Using 8i and substituting (|9|) into ([8| gives 



5i kT [Nd 
AvA = y+- + — In — 
q q \ni 
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The model uses these as the boundary conditions at ohmic 
contacts to the Si. 



B. Carrier statistics 

The equilibrium charge dens ity in the solution is given by 
the Gouy-Chapman theory. ^^^^ For simplicity, the case of a 
symmetric electrolyte will be used with the ionic species hav- 
ing a charge of either plus or minus one. Using Boltzmann 
statistics, the net charge concentration in the solution is given 

by 



C = Col 
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where Co is the salt concentration in units of cm^ . and Ef 
are assumed to be zero in the bulk of the solution. The first 
term gives the cation concentration and the second term gives 
the anion concentration. 

The equilibrium carrier statistics in non-degenerate Si are 
well known and can be stated as^^ 



n = rii exp 
p = n, exp 



( Ef-E, 
V kT 
Ej — Ef 
kT 



(13a) 
(13b) 



When modeling a semiconductor device composed of a single 
material the intrinsic energy level is identified with the elec- 
tric potential.lSHIlllin the case of heterojunctions, the discon- 
tinuity in the intrinsic levels must be taken into account.E^f^S 
In the isolated phases, the intrinsic level in the Si is £, = 5, 
relative to the Fermi (intrinsic) level of the solution, Fig.[T^). 
When the phases are brought into contact, the Si intrinsic level 
is further modified by y/, and it becomes Ej = —qy/ + 5i. If a 
voltage is applied to the Si, the Fermi level at the ohmic con- 
tact is raised by the negative of that voltage. The new Fermi 
level is Ef = —qV. Putting these values of Ei and Ef into ( 13 1 
gives 



n = rii exp 



p = n,exp 



q{y/-5i/q-V) 
kt 

q{V + Si/q-w) 
kt 



(14a) 
(14b) 



The boundary conditions for the electric potential at an ohmic 
contact can obtained assuming electrical neutrality so that, ap- 
proximately, n — Nd for n-type Sior p— Na for p-type Si. Us- 
ing this in ( [T4| gives ( [TT| , agreeing with the previous section. 

That the Fermi level is constant throughout the Si, with an 
applied voltage, can be seen from the following. In this non- 
equilibrium case, the Fermi level is replaced by quasi-Fermi 
potentials Ef —q^n for electrons and Ef — > —q^p for holes. 
The electron and hole currents are proportional to the gradient 
of the quasi-potentials.^^ Due to the insulating Si02 layer, 
the current is essentially zero, so the quasi-potentials will be 
constant across the Si and equal to each other. 



C. Diffuse layer structure and surface charges 

The Gouy-Chapman-Stern mo del wi ll be used to describe 
the diffuse layer (or double layer).™ This model states that 
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there is a thin layer, the Stern layer, on a surface which con- 
tains no salt ions from the solution. Outside of the Stern layer 
the ion concentration is given by ( [T2] l. In this paper, the thick- 
ness of the Stern layer will be set to 0.5 nm. The capacitance 
of the Stern layer is then determined by its dielectric con- 
stant. Experimental results show this capacitance to be about 
20 |jF-cm^^,'^implying a dielectric constant for the Stern layer 
of e = 1 X \Q-^^¥-cm-\ 

There can also be specifically bound charges at the Si02 
surface within the Stern layer, as discussed by Sandifer."* To 
model this a site binding model for hydroxyl groups will be in- 
cluded. Other specifically adsorbed molecules could be added 
in a similar way, for example amine groups. 

The charge on the oxide surface is due to the species MOH, 
MO^, and MOHj^, where M represents Si for a Si02 surface. 
The charge per unit area is given by ' ' 



qN, 



1 + (flf^+//:„)exp(-j3 V/) + {Kh/a^^exp{^\ir) ' 



(15) 

Ns is the site density on the oxide surface, )3 = q/kT, fl^+ is 
the activity of protons in the solution bulk, and the equilibrium 
constants are given by 



MOH+] 



and Kh = 



\M0 



(16) 



[MOH] 

The surface charge will be modeled as a uniformly charged 
thin layer with a thickness of f = 0.1 nm, at the edge of the 
Si02, with a charge density of p = a /t. 



D. Membrane 

A simple membrane model will be used to represent a layer 
of charged molecules.^' It is assumed that the charge is dis- 
tributed evenly throughout the membrane. The salt ion con- 
centration in the membrane is given by the Boltzmann distri- 
bution in the same way as in the Gouy-Chapmann model for 
the solution ( [T2] l. The charge in the membrane is then 



p =^C,„exp 



qiSm/q-yf) 



kT 



qC„,exp 



qiw-^m/q) 



kT 



(17) 

C,n is the equilibrium concentration of the salt in the mem- 
brane. It can be expressed by the partition coefficient = 
Cm/Q)^ which gives the ratio of concentration of each ion 
crossing the boundary from the solution into the membrane.^'* 
8,„ is the difference between the real potential of a solvated 
electron in the membrane and the real potential of a sol- 
vated electron in the solution, analogous to 5, in the Si. If 
the membrane is composed mainly of water, 5„, is proba- 
bly equal to zero. p„, is the uniform charge density due to 
the biomolecules. An example would be a lattice of DNA 
molecules attached to the Si02 by linker molecules as de- 
scribed in Ref.[3T] The DNA molecules have a charge of either 
1 or 2 electrons per base unit depending on whether the DNA 
strand is hybridized or not. The linker molecules are assumed 
to be uncharged and only serve to space the membrane a short 
distance away from the Si02 surface. 



IV. FINITE ELEMENT SOLUTION 

Solutions of this model require solving Poisson's equa- 
tion (|7 1, where the charge density is given for each region by 
( 12 1, ( 14 1, ( 15 I, and ( 17 1. It would be straight forward to also 



include charge in the Si02, if desired. One advantage of the 
finite element method is that the exterior boundary conditions 
are defined in a natural way and that the conditions at interior 
boundaries are matched automatically. A free, open source fi- 
nite element solver, Freefem-H-,'^ was used. In order to use 
the finite element method the Poisson equation, Qjjnust be 
converted into a variational or weak formulation P^^^ While 
this could be done in three dimensions, the examples in this 
paper will be two dimensional. 

The full length of a cylindrical nanowire can be simulated 
in cylindrical coordinates, in two dimensions using radial, r 
and axial, z, coordinates. The weak form of Q is 



/{ 



dv dw 
or or 



{vrqp} drdz — 0, 



(18) 

where v is a test function. As will be shown below, it is of- 
ten only necessary to simulate the central cross-section of a 
nanowire. In this case cylindrical symmetry is not needed and 
the weak form can be used in x-y coordinates 

^e^ + ^e^}dxdy- {vqpjdxdy^O. (19) 
ox ox oy oy ) J 

On boundaries with Dirichlet boundary conditions (ohmic 
contacts) the potential, must be specified. For the rest of 
the boundaries, Neumann conditions, the above equations as- 
sume that the perpendicular electric field is zero (by omitting 
a possible one dimensional integral on the boundary). This 
implies that there is zero current crossing the boundary. 

Since the charge density is a nonlinear function of y/, the 
above equations cannot be solved in a single step. A Newton 
iteration scheme was used^^ (see Appendix A). 

V. CALCULATIONS 

In this section, calculated results will be given in order to 
demonstrate some features of the above model. Experimen- 
tally, it is difficult to determine when the Si is at equilibrium 
with the solution. However, determination of the flat band 
condition is possible.^ The flat band voltage, Vfb, can be ob- 



tained from ( 1 1 1 by setting Ay/ = 0, for example, 

4 



q q \ni J 



p - type. 



(20) 



By expressing results relative to the flat band voltage, accurate 
values of 5,- and are not needed. 

Results shown here will use the following parameters. The 
relative dielectric constants of the solution, Si02, and Si were 
78.5, 3.9, and 11.8, respectively. The temperature was 300 K 
and the Si intrinsic carrier density was 1.45 x 10^". The Si 
is p-type with a doping concentration of 1 x 10^^ and a fixed 
hole mobility of 130cm^V^'s^'. The solution concentration 
was chosen to be 0.01 M. 
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FIG. 2: Plot of Atlas data of potential along the axis of a 0.5 [im long 
nanowire. The nanowire to solution voltage is (Vfb — 0.2)V and the 
drain-source vohage is V^s = 0.01 V. The inset shows a half cross- 
section of a 0.1 |jm long nanowire, with constant potential lines. The 
axis of the nanowire is at the bottom of the drawing. The source and 
drain contacts are at the left and right ends of the Si. The spacing of 
the potential contours is .03 V, where the line with the shortest dashes 
is -. 1 8 V and the line with the longest dashes is -.03 V. 



A. Contacts and drain current 

Contact effects and a method to calculate the device re- 
sponse can be shown with a simple cylindrical model. This 
consists of a Si nanowire of radius 10 nm with a 1 nm thick 
Si02 outer layer in a salt solution with no other charges such 
as biomolecules or site binding. There are ohmic contacts on 
both ends of the cylinder. (The Freefem-H- code is given in 
Appendix B.) 

The inset of Fig. [2] shows lines of constant potential when 
the source and drain voltages are V = (Vfl, — 0.2) V relative to 
the solution. At this potential the nanowire is in fairly strong 
depletion. The structure near the contacts does not change 
as the nanowire is made longer. The end effects due to the 
contacts extend into the nanowire a distance of roughly twice 
the radius. This end effect is independent on the length of the 
nanowire and the results were similar for other gate biases. 

While the basic response of a bioFET is the change in car- 
riers in the Si, a typical experiment will probe this by apply- 
ing a small drain-source voltage, V^s, and measuring the drain 
current, Ij. A drain-source voltage, V^s, was applied by setting 
the drain voltage (righthand contact in Fig.|2|i to V + Vds/2 and 
the source voltage to V — yds/2. When there is a current in the 
Si the above equilibrium model cannot be used. Therefore, in 
this section, current calculations were done with a commercial 
semiconductor simulation program, AtlasP^l An Atlas calcu- 
laton of the potential along the axis of the nanowire is plotted 
in Fig. |2] The applied V^s ~ 0.01 V can be seen by the dif- 
ference in y/ at either end of the nanowire. One might expect 
that the potential difference would be distributed evenly along 
the length of the nanowire with a slope of 0.02 V|jm^^. How- 
ever, the Atlas result shows that the potential curve is nearly 
flat (slope =1.67 x 10^^ V-|jm^'), except near the ends of the 
nanowire. This is because the gate (solution) essentially pins 
the potential distribution in the nanowire away from the ends. 
The current is mainly diffusion current rather than drift cur- 
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FIG. 3: (a) Plot of Atlas data of drain current vs drain-source voltage, 
at V = (Vfb — 0.2) V. The drain current is normalized by dividing by 
the value at Vjs = 0.01 V. The actual current values for Vd, = 0.01 V 
are given in the key. (b) Resistance from the Atlas data vs the length 
of the nanowire. 



rent. It would not be correct to use a gradual channel approx- 
imation with a linear voltage drop along the nanowire. 

The drain current, for three nanowire lengths, is given in 
Fig.[3];a). When the currents are normalized, the values for the 
three lengths agree. The data is linear up to about 0.02 V. The 
inverse of this slope gives the nanowire resistance, R, which is 
plotted in Fig.|3jb) for a number of lengths. It is a straight line 
with a small negative intercept. R can be accurately modeled 
as the sum of a correction due to the end resistances plus the 
resistance per unit length times the length of the nanowire. 

The central cross-section of the nanowire can be simulated, 
using the above equilibrium model and Freefem-H-. A resis- 
tance can be calculated from the average density of holes us- 
ing R = L/qptpL, where L is the length of the nanowire and 
Pt is the total number of holes integrated over the area, per 
unit length. It agrees with the resistance found by the Atlas 
simulations. This shows that it is only necessary to model the 
center cross-section of the nanowire. and that, for small Vds^ 
the full semiconductor simulation is not needed. 

One would not expect electrons to contribute to a majority 
carrier hole device with p-type ohmic contacts; the opposite 
is sometimes assumed.!^ Atlas calculations confirm that the 
electrons do not contribute, even in the case of deep deple- 
tion where the number of electrons is similar to the number of 
holes. 



B. Planar geometry 

One way to make a biofet is to use a Si on insulater 
(SOI) substrate and fabricate a ribbon shaped MOSFET.0 
A cross-section of this can be modelled. (The Freefem-H- 
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FIG. 4: Plot of parameters at the Si/solution interface of a ribbon 
geometry bioFET on a SOI substrate. The Si layer is on the left, 
with the lightly shaded band being the Si02. The hurried oxide is 
to the left of the Si and is not shown in this figure. To the right of 
this is solution, including the darker shaded DNA membrane. The 
dashed lined show the potential and hole or net ion concentrations 
at flat band bias with no DNA present. The other lines show these 
quantities at the same bias after a DNA membrane has been added. 



code for this example is given in Appendix C.) The model 
assumes a 150 nm thick hurried oxide (SOI wafer) wiht a 
20 nm thick layer of Si, doped p-type at 1 x lO'** cm^^, with 
a 1 nm oxide layer The salt solution has a concentration of 
1 X 10"2jnoles/L. 

Two situations were modelled. One with no DNA present 
and with the Si layer and the substrate bias at the same poten- 
tial to give the flat band situation. Then, with the same bias, 
a DNA layer was added. (The specific values used are given 
in the code in the appendix.) The results are plotted in Fig. |4] 
which shows only the Si layer and the region of the solution 
near the interface. With no DNA, the potential is flat accross 
the whole system and there is essentially no net charge in ei- 
ther the Si or the solution. After the DNA is added the solu- 
tion responds by having a large excess of positive ions which 
largely shield the DNA negative charge. However, there is still 
a small response in the Si, seen as the number of holes above 
the neutral 1 x 10^^ cm^^ concentration. Note that a more re- 
alistic simulation would include site binding charges on the 
oxide surface. Code for this is included in the appendix. 



C. Sensitivity 

This section will discuss the response and sensitivity of a 
nanowire to external charge as function of bias. A nanowire 
with a circular cross-section surrounded by a charged mem- 
brane was simulated. Note that there are no contacts to the Si 
on the central cross-section. Therefore, it would be wrong to 
fix the potential at the center of the nanowire as some workers 
have done.'^ The Si nanowire radius was 20 nm. The Si02 
layer was 1 nm thick and the space between this and the inner 
surface of the membrane was 1 nm. The membrane thickness 
was 3.4 nm and its fixed charge density was —4 x 10^"^ cm^^. 




0.1 0.2 0.3 
Vgate (V) 



0.5 0.6 



FIG. 5: Response, Ap, vs Vgate for a circular cross-section is plotted 
on the left axis. Sensitivity, Ap/ p(0), vs Vgate is plotted on the right 
axis. The gate voltage is given by Vgate = ^(Y ^ Vfb). 



The change in the number of holes, per unit volume, after the 
membrane is added, Ap, is plotted against gate bias in Fig.|5] 
For negative values of Vgate the nanowire is in accumulation. 
For positive bias, the nanowire is in depletion and the response 
drops off exponetially. This corresponds to the subthreshold 
regime,^ where the depletion region reaches the center of the 
nanowire. Note that this use of the term subthreshold is differ- 
ent than the common usage with respect to an inversion mode 
FET.22 

Some authors define the sensitivity as the response divided 
by the original number of holes, Ap/p{Q), which is also plot- 
ted in Fig.|5] This sensitivity increases as the depletion deep- 
ens and flattens out in the subthreshold region, in agreement 
with the approximate, analytic results of GacP and the simu- 
lation results of Liu.^ The value of this plateau depends on 
the doping level, the nanowire radius, and the concentration 
of the solution, as well as the amount of charge in the mem- 
brane. It is important to choose the best gate voltage to ob- 
tain optimal properties of an experimental device.^ Whether 
Ap or Ap/p{0) is a better indication of sensitivity will de- 
pend on the specific experimental situation. Near flat band one 
would measure a larger absolute change in current, whereas in 
the subthreshold region, the current is smaller but the relative 
change is larger 



VI. DISCUSSION 

Calculations of other effects on response can easily be done. 
Results are not shown here, but generally agree with other 
publications, for example the strong effects due to screening 
of the ions in solution^ and screening due to the site binding 
charge.^ In fact, the hole and electron densities, the ion con- 
centration, and site binding charge are strongly interdependent 
through their dependence on y/, so it is important to include 
the entire system in the same calculation. Also, if a metal- 
lic boundary is used,'^ the calculated response is much higher 
than for a semiconducting nanowire. 

Cross-sections of nanowire other than circular can be mod- 
eled. For example, calculations show that round nanowires 
and ribbons have similar sensitivity. Some other results of cal- 
culations made with this model have been published. This 
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includes the response of a trapezoidal nanowire, the effects of 
back-gating on a circular nanowire. As well as modeling pH 
measurements based on the site binding model, which agreed 
with experimental measurements. Note that the salt concen- 
tration also affects the site binding charge. 

It was also found that assuming metal boundary condi- 
tions for a nanowire gives much larger sensitivity than for the 
proper semiconductor boundary conditions. This model could 
also be used to study bioFETs in the inversion mode rather 
than the depletion/accumulation mode discussed here. In this 
case, the Boltzmann statistics for holes and electrons must be 
replaced by the Fermi distribution. One could also investigate 
partial coverage of the bioFET by charged molecules using a 
three dimensional calculation. 



VII. CONCLUSION 

A semiconductor based model of a bioFET has been de- 
veloped based on the concept of the absolute electrochemical 
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Appendix A: Newton iteration 

The finite element problem cannot be solved because the 
source term is nonlinear in i//. The approach used here is to 
start with the weak form of the problem and apply a Newton 
iteration to the integrand. 

I ^dv ^^^f ^dv ^d^f\^ 

J \dx dx dy dy J 

-j {vq[niexp{Q{V-yf))-Na 

+ Co(exp(-ev^) -exp(e\/)]} dxdy = 0. 

Where Q — q/kT, the acceptor concentration is A^^, and there 
are only holes in the Si not electrons, y/ is the potential and v 
is a test function. The first half of the second integral (second 
line) applies only to the Si and the second half (third line) 
applies only to the solution. 

Suppose the potential is given by y/ — > w + m, where wis 
a (known) guess of the potential and m is a small correction 
(to be found). Substituting yr w + m in the above equation 



gives 

/■('(Jy dw^dv 
J \dx dx dy dy J 
^ f ( dv ^du ^dv ^du^ 
J \dx dx dy dy J 

- J {vq[niexp{Q{V-w)){l-Qu)-Na+]}dxdy 

-J {[Co(exp(-ew)(l - Qu)~exp{Qw){l + Qu)]} dxdy 
= 0. 

Where exp(—Qu) was linearized using a Taylor expansion so 
that exp(— Qm) ~ (1 — Qu). Electrons can be added to the Si 
using the same method as was used for holes in the above. 
This can be rearranged to 

f ( dv du ^ dv du 1 
J \dx dx dy dy J 

- j {—qQniexp{Q{y — w))uv} dxdy 

- j {qCo [-gexp(-gH') - 2exp(2w)] uv} dxdy 

^ f ( dv dw ^ dv dw 1 
J \dx dx dy dy J 

- J {q{n,exp{Q{V - w)) ~Na) v} dxdy 

- J {qCo [exp(-gvi') - exp(gM')] v} dxdy = 0. 

The first three lines are the bilinear terms and the last three 
lines are linear Notice that —Qniexp{Q{V — w)) in the second 
line is just derivative of the hole concentration «, exp(2(V — 
w)) in the fifth line. When site binding is added, it can be 
treated the same way. The linear term will contain the site 
binding charge density, do and the bilinear term will contain 
the derivative, dao/dxj/, which gives 

- J {aov} dxdy 

Go is a complicated expression. 

Solving this finite element problem givesan approximate 
solution to M. Then a new function w + m is an improved guess 
for the true potential. Repeated iteration can further improve 
this and will tend to converge towards the true potential solu- 
tion (depending on a reasonable first guess). 
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Appendix B 

The following is an example input file to Freefem++ of a model of a cylindrical nanowire. 
// cyl02.edp 

// cylindrical nano wire MD 17 Jun 2010 

// uses Newton iteration and includes electrons and holes 
// Stern layer is commented out 
real R = 5.0e— 6; // 5e— 6 cm = 50 nm 

real Rs = l.OOOOe-6; // 1.2444e-5; //radius of the Si 

real to = l.Oe— 7; // silicon oxide thickness 

real W = l.Oe-5; 

real Q = 38.68; 

real q = 1 .6e — 19; 

real ni = 1.45el0; 

real es =1.045e-12 ;// 1 . 064 e - 12; 1 1.8*8. 85418e-14=1.04479324e-12 
real ew = 6.9505e-12; // 8.85418e-14 * 78.5,epsilon water 
//real eH = l.OOOOe-12; // 1.0000e-12= 8.85418e-14 * 11.29, Stern layer 
real esio = 3.4531e-13; // 3.9 is epsilon Si02 

real Er = —3.2528; // Si Ei minus Sol real energy , previously I used .2772 
real Na = lel8; // acceptor density , eihter Na or Nd must be zero 
real Nd = 0el7; // Donor density 

real CO = 1 .Oe -2*6.02214e20 ; // concentration in mol/L* NA / 1000 
real DO = 0el9 ; 

//real aB = l.OOOOe-7; // pH =-log_10(aB) = -In (aB)/2. 30258509299404568402 
real V = +3.5194300; // potential at ends of Si cylinder , flat band V = +3. 
real Vds = 0.0000000; // drain — source voltage , must be small . 

border Gl(t=0,Rs) {x=t; y=-W; } 
border G2( t=Rs , Rs+to ) {x=t; y=-W; } 
//border G3( t=Rs+to , Rs+to+5e-8) {x=t ; 
border G4( t=Rs+to ,R) {x=t; y=-W; } 
border G5(t=-W,W) {x=R; y=t; } 
border G6( t=R, Rs+to ) {x=t; y=^; } 
//border G7( t=Rs+to +5e -8,Rs+to ) {x=t; 
border G8( t=Rs+to , Rs) {x=t; y=W; } 
border G9(t=Rs,0) {x=t; y^; } 
border G10( t=W,-W) {x=0; y=t ; } 

border Tl(t=\Y-W) {x=Rs; y=t ; } 
border T2(t=-W,W) {x=Rs+to ; y=t ; } 
//border T3(t=-W,W) {x=Rs+to+5e-8; y=t ; } 

mesh Th = buildmesh (Gl(40) + G2(5) + G4(20) + G5(20) + 

G6(20) + G8(5) + G9(40) + G10(200) + 

Tl(200) + T2(200)); 
//mesh Th = buildmesh (Gl (20) + G2(5) + G3 ( 1 ) + G4(2) + G5 (20) + G6(20) + 
// G7(20) + G8(2) + G9(1) + G10(5) + G11(20) + G12(200) + 

// Tl(200) + T2(200) + T3(200) + T4(200)); 

fespace Ph(Th,PO); 
fespace Vh(Th,Pl); 

Ph reg = region ; 

plot(reg , fill=l,wait=l,value=l); 

int nSi=reg(le-8,0); 
int no=reg (Rs+to /2 ,0); 

//int nH=reg (Rs+to +2.5e -8 ,0); // Inner Helholtz or Stern layer 



y=-W; } // Stern layer 
y^; } 
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int ns = reg (Rs+to+2e — 7 ,0); 

cout <<", nSi "«nSi<<" ,no "« no <<", ns "«ns« endl ; 
Vh u, V, w, u2 ; 
real cpu=clock ( ) ; 

// Stern layer + eH/q*( region==nH) 

Ph epq = es /q *( region==nSi ) + esio /q * ( region==no ) + 
ew/q*(region==ns ); 

// The first guess for potenial . Must contain the correct boudary conditions. 
w = ( (V+ Vds*yAV) + Er -(1/Q)* log (Na/ ni ) ) * ( sin ( pi *y /(2 *W) ) ) ^2 *(x<Rs); // for p-type 
u = 0; 

cout « w(le-6,W) « " "« w(le-6,-W) «" "«w(0,0) «" "« w(R,0) « endl; 

problem Lap (u , v , sol ver=LU) = 

i nt2d (Th ) ( ( dx (u)* epq*dx ( V ) + dy ( u) * epq*dy ( v ) ) * x ) // Laplace 

+ int2d(Th)( ni*Q*( exp (Q*(V+Er-w))+ exp (Q*(w^Er-V)) ) * u*v*x *( region==nSi ) ) // Si 

— int2d(Th)( C0*(— Q* exp(— Q*w)— Q* exp (Q*w)) * u*v*x *( region==ns ) ) // sol 
+ int2d(Th)( (dx(w)* epq*dx( v) + dy (w) * epq*dy ( v ) ) * x ) // Laplace 

-int2d(Th)( (ni*( exp (Q* (V+Er-w)) - exp (Q* (w^Er-V) ) )-Na+Nd) * v*x * ( region==nSi ) ) // 

— int2d(Th)( CO *( exp(— Q*w)— exp (Q*w) ) * v*x * ( region==ns ) ) // sol 
+ on(G5,u=0) + on(Gl,u=0) + on(G9,u=0); 

Lap ; 

cout «"w+u, w, u: "« w(Rs,0) + u(Rs ,0) «" "« 
w(Rs+le-7,0)<<" "« u(Rs,0) « endl; 

//plot (Th, w, value = true , coef = l, wait = true); 
plot (Th, u, value = true , coef = l, wait = true); 

u2 = w + u/2; 

w = u2 ; 

// Iterate until solution converges; ie . u < about e — 16 everywhere 
for (int j=l;j<=20;j++){ 
Lap ; 

cout « j <<") w+u ,p ,u: "« w(Rs,0) + u(Rs ,0) «" "« 
ni * exp (Q* ( V+Er-w( Rs , ) ) ) < < " " «u ( Rs , 0) < < 
endl ; 
u2 = w + u/ 1 ; 
w = u2 ; 
} 

// calculate the total number of holes (per cm length) on the cross —section 

// at the middle of the cylinder 

int Nn = 80; 

real X = 0; 

real xdel = Rs/Nn; 

real pdel = 0; 

real ptot = 0; 

for (int n =1 ;n<=Nn; n ++){ 

X = n* xdel —xdel /2 ; 

pdel = 2* pi * ni *exp (Q*(V+Er— w(X, 0) ) ) * ( n*xdel— xdel /2) * xdel ; 

ptot = ptot + pdel ; 

} 

cout «"total holes "« ptot <<" pdel "« pdel «" p_ave "« ptot /( pi *Rs'^2) « endl; 
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// now calculate the total number of electrons 
X=0; 

real ndel = 0; 

real ntot = 0; 

for (int n =1 ;n<=Nn; n++){ 

X = n*xdel— xdel /2 ; 

ndel = 2* pi * ni * exp (Q* (w(X,0) — Er— V) ) * ( n* xdel— xdel /2) * xdel ; 
ntot = ntot + ndel ; 

} 

cout «"total electrons "« ntot «" n_ave "« ntot /( pi *Rs'^2) « endl ; 

Vh Ex = -dx(w); 
//func Ex = — dx(w); 

Ph unite = 1 *( region==ns ) ; 

//Ph Cout = C0*( exp(— Q*w) — exp (Q*w)) * ( region==ns ) ; // This gives digital steps 
Ph unitpNa = ( region==nSi ) ; // region of the Si trapezoid 
//Ph unitf = ( region==nf ) ; 
X = 0; 

real Ni = 400; 
{ 

ofstream f ile (" cy2 . dat " ) ; 

file « "# ave hole density over area of the Si at y = 0, "« 
ptot/( pi*Rs'^2) « "\n"; 

file « "# ave electron density over area of the Si at y = 0, "« 
ntot/(pi*Rs^2) « "\n"; 
for (int i =0; i <=Ni ; i ++) { 
X = + i *(0+R)/Ni; 

file « X « " " « w(X,0) « " " « 

unitC(X,0)*CO*(exp(-Q*w(X,0))-exp(Q*w(X,0))) «" "« 

unitpNa(X,0)*( ni *exp (Q* (w(X,0) - Er-V))) < <" "« 

unitpNa(X,0)*( ni *exp (Q*(V+Er-w(X, 0) ) ) ) « 

" "« unitpNa(X,0)*(-Na+Nd) « 

" "« Ex(X,0) « "\n"; 

} 

}; 

real Y; 
{ 

ofstream f ile (" cyz2 . dat " ) ; 
for (int i=0;i<=Ni; i++){ 

Y = -W + i*2*W/Ni; 

file « Y «" "« w(0,Y)<<" " « w(Rs/2,Y) « "\n"; 
} 

} 



/* 

Vh Ey = -dy(w); 
real Nc = 100; 
{ 

ofstream f ile (" nwc2b . dat " ) ; 
for (int i =0; i <=Nc ; i ++){ 
X = + i*Rs/Nc; 

file « X «" "« Ey(X,W)<<" " « Ey(X,-W)«" " « 
q*100*ni*exp(Q*(V+Er-w(X,0)))*Ey(X,W) « "\n"; 

} 

} 
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*/ 

plot (Th, w, value = true , coef = l, wait=true); 
plot (Th, u, value = true , coef = l, wait = true); 
u2 = w+u ; 

real[int] viso(6); 

// viso [0]=.03; 

// viso [1]=0; 

// viso [2]= -.03; 

// viso [3]= -.06; 

//viso[4]=-.09; 

for ( int i =0; i <6; i ++) 

viso [i]=i*-0.03 -0.03; 
//for (int i=4;i<9;i++) 
// viso[i]=i*-0. 020-0. 01; 
//for (int i =9; i <viso . n ; i ++) 
// viso [i]=i * -0.100+0.63; 

plot ( u2 , viso = viso (0 : viso . n — 1) , value = true , ps = "cy— pot . eps " , coef = l, wait=true); 



Appendix C 

The following is an example input file to Freefem++ of a model of a ribbon nanowire os a SOI substrate. 

// file "plane -16. edp" MD 20 Feb 2012 

// from planar goemetry with site binding and DMA 

// 

// Copyright National Research Council of Canada, 2012. 
// This is an input file for Freefem++ 
// with hydroxyl and amine site binding 
// 

real R = 1.2e— 5; // Right hand edge of the electrolyte , 5e— 6 cm = 50 nm 

real L = 1.5e— 5; // Left edge of the hurried oxide (BOX) 

real W= 1.4e— 5; // half height of the modellel region 

real ts = 2e— 6; // thickness of Si layer 

real to = le— 7; // thicknes of Si02 

real tc = le — 8; // thickness of site binding layer 

real tH = 4e — 8; // Stern layer thickness is the sum of tc+tH, ie 5 Anstroms 
real tL = 5.0e— 8; // thickness of linker beyond tH , must be => le— 8 
real tD = 3.40e— 7; // thickness of membrane, 3.4e— 7 is 10 DNA bases 
real Q = 38.68; // Q = q/kT 
real q = 1.6e — 19; // unit charge 

real ni = 1.45el0; // Si intrinsic concentration 

real es = 1.045e-12; // dielectric constant of Si, 1 1 .8*8.85418e-14 

real ew = 6.9505e-12; // 8.85418e-14 * 78.50 , epsilon water 

real eH = l.OOOOe-12; // 1.0000e-12= 8.85418e-14 * 11.29, Stern layer 

real esio = 3.4531e — 13; // 3.9 is epsilon Si02 

real em = ew; // 1 .7708 e — 12; // membrane dielectric constant 

real Na = lel8; // acceptor density , eihter Na or Nd must be zero 

real Nd = 0el7; // Donor density. Nd aded 10 Jun 2010 

real CO = 1 .Oe -2*6.02214e20 ; // concentration in mol/L* NA / 1000 

real Er = —3.7000; // energy difference Si — electrolyte (I used 0.2772 previously) 
real Cm = 1.000*C0; // ion concentration in membrane 

real Em = 0.000; // intrinsic energy level of membrane wrt electrolyte 
real cDn = 2; // cDn = 1 for ssDNA, cDn = 2 for dsDNA 

real cD = — 4.000e20*cDn ; //volume charge density of DNA strip , ie . membrane 
real Ss ; // surface density of site binding charges, calculated below 
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real Ns = ; // 1 . 5 el4/ 1 e -8; // site binding site density in 0.1 nm thick layer (5 .0 el4/l e-8) 
real Ka = le+2; 
real Kb = 1 e — 6 ; 

real Nam = 0;//l .5 el4/le— 8; // site binding site density due to amine sites 
real Kam = le — 10; // amine disociation constant 
real aB = le— 7; // pH =7 gives aB = le— 7 
// Vfb = -Er + (l/Q)*ln(Na/ni) 

real V = 4.1666266077 — 0.00; // voltage at the drain and source contacts 
// for flat band Vback = -(Er + Bwf) 

real Vback = +4.2272 — 0.00; // potential at x= — L, assumes the work function of the metal 
real Bwf = —0.5272; // Back Si wafer work function = —0.5272 for degenerate p— type 
/* Vback is applied to an ohmic contact to Si on the letf edge of the BOX 

* The boundary value of the potential must include the workfunction of 

* the Si wafer. If it is degenerate p— type this equals the valance band energy, 

* which is —5.25 eV or —0.5272 w.r.t. the Si intrinsic level. The reference level 

* in this model is the work function of the electrolyte and the Si intrinsic level 

* is Er w.r.t. that. So the left boundary condition for potential is 

* Vback + Er — 0.5272 or Vback + Er + Bwf which is used in the first guess for w. 
*/ 
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=-L,0) {x=t; y=-W; } 

=0,ts) {x=t; y=-W; } 

=ts,ts+to) {x=t; y=-W; } 

=ts+to , ts+to+tc ) {x=t; y=-W; } 

= ts+to+tc , ts+to + tc+tH) {x=t; y=-W; } 

t = ts + to + tc+tH , ts+to + tc+tH+tL) {x=t; y=-W; } 

t = ts + to + tc+tH+tL , ts + to + tc+tH+tL+tD) {x=t; y=-W; } 

= ts + to + tc+tH+tL+tD , ts + to + tc+tH+tL+tD+R) {x=t; y=-W; } 

=-W,W) {x=ts+to + tc+tH+tL+tD+R; y=t ; } 

= ts + to + tc+tH+tL+tD+R, ts + to + tc+tH+tL+tD) {x=t; y=W; } 

t=ts + to + tc+tH+tL+tD , ts + to + tc+tH+tL) {x=t; y=W; } 

t = ts+to + tc+tH+tL , ts+to + tc+tH) {x=t; y=W; } 

= ts + to + tc+tH , ts + to + tc ) {x=t; y=W; } 

t=ts+to+tc , ts+to ) {x=t; y=W; } 

t=ts+to , ts ) {x=t ; y=W; } 

t=ts ,0) {x=t ; y=W; } 

t=0,-L) {x=t; y=W; } 

t=W,-Wi) {x=-L; y=t; } 

=-W,W) {x = 0; y=t; } 

=W,-W) {x=ts ; y=t ; } 

=-W,W) {x= ts + to; y=t; } 

=W^-W) {x=ts+to + tc; y=t; } 

=-W,W) {x=ts + to + tc+tH; y=t; } 

=W,-W) {x=ts + to + tc+tH+tL ; y=t; } 

=-W,W) {x=ts+to + tc+tH+tL+tD ; y=t ; } 



mesh Th = buildmesh (Gl (40) + G2(20) + G3 (4) + G4(4) + G5 (3) + G5b(3 ) + G6a(3) + G6(20)+G7 (40) 

+G8(20) + G8b(3) + G9a(3) + G9(3) + G10(4) + Gll(4) + G12(20) + G13(40) + G14(40) 
+T1(80) + T2(80) + T3(80) + T4(80) + T5(80) + T6(80) + T7(8 0)); 
//mesh Th = buildmesh (Gl (40) + G2(20) + G3 (4) + G4(2) + G5 (3) + G6(20) + G7 (40) 
// +G8(20) + G9(3) + G10(2) + G11(4) + G12(20) + G13(40) + G14(40) + 

// Tl(80) + T2(80) + T3(90) + T4(90) + T5(9 0)); 



fespace Ph(Th,P0); 
fespace Vh(Th,Pl); 
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Ph reg = region ; 

plot(reg , fill=l,wait = l,value = l); 

int nf = reg(— L/2,0); // buried oxide 
int nt = reg ( ts /2 ,0) ; // Si 

int no=reg ( t s + to /2 , 0) ; // thin Si02 on surface of Si 

int nc=reg ( ts + to + tc /2 ,0); // charged site binding layer 

int nH=reg ( ts+to + tc+tH/2 ,0) ; // Inner Helholtz or Stern layer 

int nL=reg ( t s + to + 1 c+tH+tL /2 , 0) ; // solution in linker region 

int nm=reg ( ts + to + tc+tH+tL+tD/2 ,0); // DNA membrane 

int ns=reg (R— le — 7 ,0); 

cout «"nf "<<nf<<", nt "<<nt<<", no "<<no<<" ,nc "<<nc« 

,nH "«nH<<", nL "<<nL<<", nm "<< nm <<", ns "«ns« endl ; 

Vh u , V , w ; 

// The function ep gives the dielectric constant for each region 
Ph ep = esio *( region==nf ) + es *( region==nt ) + esio *( region==no) + 
eH* ( region==nc ) + eH* ( region==nH) + ew*( region==nL) + 
em* ( region==nm) + ew*( region==ns ) ; 

// This is the first guess for the potential , w. It must satisfy the boundary conditions. 
// That is the potential at the left edge of the BOX must equall Vback 
// and the potential at the right edge of the electrolyte must be 0. 
w = ( Vback+Er+Bwf)*(x/(-L))*(x<0); 

//u is the correction (function) in the Newton iteration. 

// This problem solves for u, which should get small, about le — 17 on convergence, 
problem Lap (u , v , s ol ver=LU) = 

int2d (Th) ( ( dx (u)* ep*dx ( v) + dy (u) * ep *dy ( v ) ) ) // Laplace 

+ int2d(Th)( q*ni*Q*( exp (Q*(V+Er— w)) + exp (Q* (w^Er— V) ) ) * u*v *( region==nt ) ) // Si 

— int2d(Th)( q*CO*(— Q*exp(— Q*w)— Q*exp (Q*w))* u*v*( region==ns ) ) // sol 

— int2d(Th)( q*CO*(— Q* exp(— Q*w)— Q*exp (Q*w)) * u* v *( region==nL) ) // sol in linker 

— int2d (Th)( u*v*q*Ns*( -(aB /Ka) *Q*exp(-Q*w) -(Kb/aB )*Q* exp (Q*w) -4*(Kb/Ka) *Q)/ 

( 1 + ( aB /Ka) * exp(— Q*w) + (Kb/aB ) * exp (Q*w))'^2 * ( region==nc ) ) // hydroxyl site 

— int2d(Th)( u*v*q*Nam*( — (aB/Kam)*Q*exp(— Q*w) )/ 

(l + (aB/Kam)* exp(— Q*w))'^2 *(region==nc) ) // amine site 
+ int2d(Th)( q*Cm*Q*( exp (Q*(Em-w))+ exp (Q*(w— Em)) ) * u* v * ( region==nm) ) // membrane 
+ int2d(Th)( dx(w)* ep*dx(v) + dy (w) * ep *dy ( v ) ) // Laplace 

-int2d(Th)( q*(ni*( exp (Q* (V+Er-w)) - exp (Q* (w^Er-V) ) )-Na+Nd) * v * ( region==nt ) ) // Si 

— int2d(Th)( q*CO* ( exp(— Q*w) — exp (Q*w) ) * v * ( region==ns ) ) // sol 

— int2d(Th)( q*CO* ( exp(— Q*w) — exp (Q*w) ) * v * ( region==nL) ) // sol in linker 

— int2d(Th)( v*q*Ns *( aB /Ka*exp(—Q*w)— Kb/aB* exp (Q*w) ) / 

(l + aB/Ka*exp(— Q*w) + Kb/aB*exp(Q*w)) *( region==nc ) ) // hydroxyl site 

— int2d(Th)( v*q*Nam* (aB /Kam* exp(— Q*w) ) / 

(l + aB/Kam*exp(— Q*w)) * ( region==nc ) ) // amine site 

— int2d(Th)( q*(Cm*( exp (Q* (Em-w)) — exp (Q* (w^Em) ) ) + cD) * v * ( region==nm) ) // Membrane 
+ on(G7,u=0) + on(G14,u=0); 

Lap ; 

cout «"w+u, w, u: "« w(ts ,0) + u(ts ,0) «" "« 
w(ts,0)<<" "« u(ts,0) « endl; 

plot (Th, u, value = true , coef = l, wait = true); 

w = w + u ; 

real Xss = ts+to+tc/2; // center of site binding charge at y = 
real Xsur = ts ; // surface of the Si at y = 
// TODO: write a test for convergence for the following loops 
// Right now I manually adjust the number of loops 
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for (int j =l;j <=ll;j ++){ 
Lap ; 

Ss = Ns*(aB/Ka*exp(-Q*w(Xss,0))-Kb/aB*exp(Q*w(Xss ,0)))/ 

(l + aB/Ka*exp(-Q*w(Xss ,0)) + Kb/aB*exp(Q*w(Xss ,0))) 
+ Nam*(aB /Kam*exp(— Q*w(Xss , ) ) ) / 
( 1 + aB /Kam*exp(— Q*w(Xss , ) ) ) ; // sum of hydroxyl and amine charges 
cout « j <<") w+u ,p ,u: "« w(Xsur ,0) + u(Xsur ,0) «" "« 
ni*exp(Q*(Er-w(Xsur ,0))) < <" "« 
u(Xsur,0)<<" sigma "« Ss*le-8 « end! ; 
w = w + u ; 
} 

// The above iterations take too long to converge, for some values of V and Vback . 

// At this point use the solution , w, to generate a new mesh 

// with highest grid density in regions where w changes most. 

// error= 0.002 controls the grid density and was found by trial and error. 

// hmax = 2.5e— 7 sets the largest grid size. This is needed to get an accurate enough 

// solution in the electrolyte region. Again, the value was found by trial and error. 

real error =0.0020; 

Th=adaptmesh (Th , w, err = error ,hmax = 2.5e— 7,nbvx=38000); 
for (int j=l;j <=28;j++){ 
Lap ; 

w = w + u ; 

Ss = Ns*(aB/Ka*exp(-Q*w(Xss,0))-Kb/aB*exp(Q*w(Xss ,0)))/ 

(l + aB/Ka*exp(-Q*w(Xss ,0)) + Kb/aB*exp(Q*w(Xss ,0))) 
+ Nam*(aB/Kam*exp(— Q*w(Xss , ) ) ) / 
( 1 + aB /Kam*exp(— Q*w(Xss , ) ) ) ; // sum of hydroxyl and amine charges 
cout « j <<") w ,p ,u: "« w(Xsur,0) «" "« 
ni*exp(Q*(V+Er-w(Xsur ,0))) < <" "« 
u(Xsur,0)<<" sigma "« Ss*le-8 « endl ; 

} 



Vh Ex = -dx(w); 

//func Ex = — dx(w); 

Ph unite = 1 *( region==ns ) ; 

Ph unitL = 1 *(region==nL) ; 

Ph unitm = ( region==nm) ; // region of the membrane 

//Ph Cout = C0*( exp(— Q*w) — exp (Q*w)) * ( region==ns ) ; // This gives digital steps 
Ph unitpNa = ( region==nt ) ; // region of the Si trapezoid 
//Ph Dout = w*( region==nt ) ; 

// calculate number of total holes , using Freefem++ built in integration 

// This is for a 1 cm long trapezoid (nanowire). 

real Ptot = int2d(Th)( ni * exp (Q* (V+Er— w) ) * ( region== nt ) ); 

cout «"total holes integrated over area of Si "« Ptot « endl; 

cout «"average holes "«Ptot /( ts *2*W)<< endl; 

// calculate number of total electrons , using Freefem++ built in integration 

real Ntot = int2d(Th)( ni * exp (Q* (w^V— Er ) ) * ( region==nt ) ); 

cout «"total electrons integrated over area of Si "« Ntot « endl; 

cout «"average electrons "<<Ntot /( ts *2*W)<< endl; 

// The area of the Si is Area = ts *(2*W) 
real Areaint = int2d(Th)( (region ==nt) ); 

cout «"Area "« ts *(2*W)<<" , Arealnnt "«Areaint « endl; 
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cout <<"total surface charges (per cm'^2) "<< Ss * le— 8 « endl ; 

cout «"hydroxyl charges (per cm'^2) "<< Ns *( aB /Ka* exp(-Q*w( Xss ,0)) - Kb/aB*exp (Q*w(Xss , ) ) ) / 
(l + aB/Ka*exp(-Q*w(Xss ,0)) + Kb/aB*exp(Q*w(Xss ,0))) * le-8« 
" , amine charges (per cm'^2) "«Nam* (aB/Kam* exp(— Q*w(Xss , ) ) ) / 
(l + aB/Kam*exp(-Q*w(Xss ,0))) * le-8 « endl; 

real X; 

real Ni = 1350; // number of data points in the file "pi— 16.dat" 

{ 

ofstream f i le (" pi - 16. dat " ) ; 

file «"# Vback = "« Vback <<", V = "« V 

<<", Area of Si (cm'^2)= "« ts*(2*ts) « "\n"; 
file <<"# total holes in Si layer of 1 cm length = "<< Ptot « "\n"; 
file <<"# total electrons in Si layer of 1 cm length = "<< Ntot « "\n"; 
file « "# the surface charge (at Xss="« Xss <<") is (l/cm'^2) 
" « Ss*le-8 « "\n"; 

file « "# Data for the following parameters is given along Y=0."<< "\n"; 
file <<"# X potential electrolyte membrane electrons holes dopant 
electric_field "<<"\n"; 
for (int i =0; i <=Ni ; i ++){ 
X = -L + i *(L+ts + to + tc+tH+tL+R)/Ni ; 
file « X « " " « w(X,0) « " " « 

unite (X,0) * CO*(exp(-Q*w(X,0)) - exp (Q*w(X, ) ) ) 

+ unitL(X,0)*CO*(exp(-Q*w(X,0))-exp(Q*w(X,0))) «" "« 

unitm(X,0)*( Cm*( -exp (Q* (w(X,0) -Em) ) + exp (Q* (Em-w(X, ) ) ) ))<<" "« 

unitpNa(X,0)*( ni *exp (Q*(w(X,0) -Er-V))) < <" "« // electron density 

unitpNa(X,0)*( ni *exp (Q*(V+Er-w(X, 0) ) ) ) « 

" "« unitpNa(X,0)*(-Na+Nd) « // Nd added 10 Jun 2010 

" "« Ex(X,0) « "\n"; 

} 

}; 

real N; 
real Y; 

real Tspace = 2e— 7; // spacing of data points , 2e— 7 = 2 nm 
{ 

ofstream f i le (" tr —ex— pot . dat " ) ; 

file «"# Potential along Si surface"« "\n"; 

file «"# Vback = "« Vback <<", V = "« V <<", site density = "« Ns « "\n"; 
file <<"\n\n# X, Y, and Potential along the front surface of the Si"« "\n"; 
N = (W)/Tspace ; 
N = 2*round(N); 

for (int i=0;i<=N; i++){ 
X = ts ; 

Y = (-W) + i *2*(W)/N; 

file « X «" "« Y «" "« w(X,Y) « "\n"; 

} 

}; 

plot (Th, u, value = true , coef = l, wait = true); 
plot (Th, w, value = true , coef = l, wait=true); 
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